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ABSTRACT 

NIMASTEP is a dedicated numerical software developed by us, which allows one to integrate the osculating motion 
(using cartesian coordinates) in a Newtonian approach of an object considered as a point-mass orbiting a homogeneous 
central body that rotates with a constant rate around its axis of smallest inertia. The code can be applied to objects 
such as particles, artificial or natural satellites or space debris. The central body can be either any terrestrial planet 
of the solar system, any dwarf-planet, or even an asteroid. In addition, very many perturbations can be taken into 
account, such as tiie combined third-body attraction of the Sun, the Moon, or the planets, the direct solar radiation 
pressure (with the central body shadow), the non-homogeneous gravitational field caused by the non-sphericity of the 
central body, and even some thrust forces. The sinnilations were performed using different integration algorithms. Two 
additional tools were integrated in the software package; the indicator of chaos MEGNO and the frequency analysis 
NAFF. NIMASTEP is designed in a flexible modular style and allows one to (de)select very many options without 
compromising the performance. It also allows one to easily add other possibilities of use. The code has been validated 
through several tests such as comparisons with numerical integrations made with other softwares or with semi-analytical 
and analytical studies. 

The various possibilities of NIMASTEP are described and explained and some tests of astrophysical interest are 
presented. At present, the code is proprietary but it will be released for use by the community in the near future. 
Information for contacting its authors and (in the near future) for obtaining the software are available on the web site 
http : //www . f undp . ac . be/en/research/pro j ects/page_view/ 10278201/ 

Key words. Methods: numerical - Gravitation - Celestial mechanics - Space vehicles 



1. Introduction 

The orbital motion of an artificial satellite, a space debris, or 
a natural satellite have important similarities: their dynam- 
ics mainly depends on the physical parameters of the cen- 
tral body (mass, radius, shape, etc.), the physical param- 
eters of the object under study (mass, shape, reflectivity, 
etc.) and the dynamics of other bodies orbiting the same 
central body (for example the influence of the Moon or 
of the Sun on an Earth's satellite). Therefore, if the ini- 
tial orbital conditions of the considered body (semi-major 
axis, eccentricity, inclination, etc.) and the parameters men- 
tioned above are known, the computation of the orbital dy- 
namics can be performed, at a certain level of accuracy, 
using the same program whatever the exact nature of the 
problem. 

Of course, the orders of magnitude of the forces act- 
ing on the body will difl'er from case to case. For example, 
the direct solar radiation pressure will have a significant 
effect on an artificial satellite of the Earth but none on a 
natural satellite of an asteroid. In contrast, the effect of 
the non-sphericity of the central body will have a much 
greater influence on a satellite of asteroid (with a shape 
usually far from a sphere) than on an artificial satellite of 
the Earth. Therefore an incomplete dynamical model with 
forces switched on/off' depending to magnitude can be use- 
ful to explain the main part of the real motion or to validate 



an analytical theory (in which we keep only the main ef- 
fects). 

These considerations have led us to create the soft- 
ware NIMASTEP (Numerical Integration of the Motion 
of Artificial Satellites orbiting a TElluric Planet). 
NIMASTEP has been (first) implemented during the Ph.D. 
thesis of N. Delsate (Delsate 2011b). NIMASTEP was ini- 
tially developed to numerically integrate the motion of an 
artificial satellite orbiting a telluric planet (Delsate et al. 
2010; Lemaitre et al. 2009; Valk et al. 2009a)i. Later on this 
software was extended to calculate the motion of one or sev- 
eral artifical (Delsate 2011a) or natural satellites (Compere 
et al. 2012) orbiting a (dwarf-)planet or an asteroid of the 
solar system. In the current version of the software, several 
forces are still neglected. This will be discussed briefiy in 
section 2.3. But the current state of the code is not the 
limit, future versions could include these forces step by 
step. The aim of this software is not to become an orbit 
restitution software but to be used to analyze the different 
orbital dynamics. So far, the code can handle quite a few 
different applications. The first step is always a two-body 
problem, where the central body is considered as a homo- 
geneous mass and rotates (at constant rate) around its axis 
of smallest inertia. The small body is distant and small 
enough to have no effect on the motion of the central body 



^ The name of this software has been fixed at this time. 
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and spherical enough to be considered a point-mass satel- 
lite. As we will show below, different pertubations already 
available can be applied to this initial system. 

The first characteristic of our software is its extensibility 
and its generality (within its possibilities) : it allows a huge 
choice of central bodies, of forces, of objects of interest, 
etc. and can therefore be used by various members of the 
dynamical or geodesian community. The second originality 
of our software is the availability of two efficient tools of 
analysis: the indicator of chaos MEGNO (Cincotta et al. 
2003; Cincotta & Simo 2000; Gozdziewski et al. 2001) and 
the frequency analysis NAFF (Laskar 1999, 2003, 2005) are 
directly connected to the numerical integrations. 

The code is written in FORTRAN 90 and is divided 
into modules that can be used in other programs (e.g., the 
module of the numerical integrators). Each force (gravity 
field of the central body, radiation pressure, gravitational 
attraction of a third body, etc.) can be switched on or off. 
The software has been validated through several methods: 
comparison with other integrators, checkings with analyt- 
ical and semi-analytical existing methods and systematic 
reproduction of known behaviors. 

The paper is structured as follows: in section 2, we de- 
scribe the forces that can be taken into account (or not) 
by NIMASTEP and we give the related equations of mo- 
tion. In section 3, we present the two supplementary tools 
(MEGNO and NAFF). In section 4, we demonstrate the re- 
liability of NIMASTEP by reproducing known results; we 
present some tests to show that NIMASTEP correctly sim- 
ulates physical processes by comparisons with another soft- 
ware and with semi-analytical results and analytical stud- 
ies. In the last section, we summarize and highlight the 
major features of the software and its future extensions. 

The version of NIMASTEP presented in this paper is 
the web version 5.35w. 



2. Code capabilities 

NIMASTEP numerically integrates the osculating equa- 
tions of motion in a Newtonian approach of a small body 
considered as a point-mass around a homogeneous and ro- 
tating (with a constant rate) main body. The (fixed) refer- 
ence frame has its origin at the center of mass of the main 
body and the x-y plane corresponds to the equatorial plane. 
The location of the prime meridian (sidereal angle) and the 
direction of the north pole (right ascension and declination) 
with respect to the International Celestial Reference Frame 
(ICRF) are extracted from Archinal et al. (2011). 

The position of the small body is written in carte- 
sian coordinates {x,y, z,x,y, z)={r,r = v) and r = \r\ = 
y/x^ + y'^ + z^. This choice ensures a high stability for the 
numerical solution by avoiding both eccentricity and incli- 
nation singularities. The equations of motion are given by 

' r = V 

< . .. ^ (1) 

V = rpot + 2^ T-36, + vvrp + rth , 

^ i 

where rpot, i^'ibi, vr-r-p and rth are the accelerations due to 
the gravitational potential of the central body, respectively 
to the gravitational attrac'tion of the third bodies, to the 
direct solar radiation pressure (with or without umbra) and 



to the possible thrust acting on the satellite. These accel- 
erations are developed in section 2.2. 

The integrations are performed in the fixed reference 
frame (defined above). However, the computation of some 
parts of the acceleration, as rpot, is made in a uniformly 
rotating frame (of velocity 6 around the axis of smallest 
inertia of the main body). To convert it to the fixed frame, 
we use 

rflx = Re r^ot , (2) 

with rfix (rrot) the location of the satellite in the fixed refer- 
ence frame (in the rotating reference frame), 9 the sidereal 
angle (time-dependent) and Re the rotation matrix given 
by 

/cos6l -sin6l 0\ 
Rg = sin6l cose . (3) 
V 1/ 

Equation (2) corresponds to the following expression for 
the acceleration: 

rfix = Re rrot + 2 -Rg r^ot + Re r^t , (4) 

where the values of 9 and 9 come from Archinal et al. 
(2011). We consider 9 to be equal to 0. 

2.1. Numerical integration 

The software enables the user to choose the integrator. 
The present propositions are on the one hand three single- 
step Runge-Kutta methods (hereafter RK): a Runge-Kutta 
of order 4 (Hairer et al. 1993) (RK4), a Runge-Kutta- 
Fehlberg of order 4 with a variable step size (Hairer et al. 
1993) (RKF4), and an explicit Runge-Kutta method of 
order 8(5,3) (DOP853) according to Dormand & Prince 
(Hairer et al. 1993); on the other hand, some other methods: 
a Bulirsch-Stoer (Stoer & Bulirsch 1980) with a variable 
step size, and the multistep method of Adams-Bashforth- 
Moulton of order 10 (Hairer et al. 1993) (ABMIO). Each of 
these integrators has been tested and validated. 

The first method (RK4) is one of the most widely used 
methods. It is efficient for simple problems with smooth 
curves. But it is not adapted to problems with high vari- 
ation of the derivatives, such as highly eccentric orbits. 
The RKF4, DOP853, and Bulirsch-Stoer use an adapta- 
tive step size to improve, for example, the integrations of 
these orbits. Runge-Kutta method of order 8 (DOP853) and 
Bulirsch-Stoer are algorithms often used in spatial geodesy 
and n-body problem integrations. This is why we intro- 
duced these integrators in our software. Finally the ABMIO 
was implemented to obtain a precise constant-step-size in- 
tegrator faster than RK methods. 

The advantage of multistep over single-step RK meth- 
ods of the same accuracy is that the multistep methods 
require only one function evaluation per step, while, e.g., 
the RK4 method requires four, and the RKF4 method, six 
function evaluations. Conversely, the disadvantage of the 
multistep methods is the difficulty to change the step size 
during the integration, while for the single-step RK meth- 
ods this is an easier procedure. The main disadvantage of 
explicit methods (like RK) is that they are only condition- 
ally stable, and the critical time step may be several or- 
ders of magnitude smaller that the desired time step. The 
ABMIO overcomes this disadvantage by using first an ex- 
plicit predictor and second an implicit corrector (uncondi- 
tionally stable). The advantage of Bulirsch-Stoer is that it 
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is often possible to choose extremely long step sizes while 
maintaining reasonable accuracy. But the algorithm is more 
complicated and makes a single step rather slow. 

2.2. Accelerations 

Let us describe the different accelerations (or forces) taken 
into account. 

2.2.1. Non-sphericity of the central body: spherical 
harmonics expansion 

The calculation of the acceleration through the gravity field 
of the central body is performed by using a spherical har- 
monics development for the potential (Kaula 1966) 



oo n / T-) \ n 



C/(r,A,0) = -^^^ 



r ^ — ' ^ — ' V r 

n=0 m=0 



P„,„(sin^) X (5) 



(Cn,m cos(mA) + Sn^m sin(mA)j , 

where n = GM is the gravitational constant of the central 
body, -Re is the radius of the main body, (r, A, 0) are the 
spherical coordinates of the small body, Cn,m and Sn^m are 
physical constants depending on the shape of the main body 
and named coefficients of the expansion (n is the degree 
and m is the order of the coefficients) and Pn,m are the 
associated Legendre polynomials. 

To calculate the partial derivatives of the potential U = 
H V with respect to the cartesian coordinates, and then the 
acceleration, we use the accurate and efficient formulation 
(in the rotating frame) of Cunningham (1970): 



.. _ (9V_ dV_ dV_y 
where ^ = 3? < I] I] i?e (^ 



(6) 



n=0 m=0 



dm 



with 5R {z} giving the real part of the complex z. The inter- 
mediate functions 14,,™ are calculated with the recurrent 
formulas 



^0,0 


= 


(7) 




= (2n - 1)^:2 — K-l,n-l , 


(8) 


Vn+l,n 


= (2n-hl)^ F„,„, 


(9) 


(n - m)Vn,m 


= (2n 1) K-i,™ 


-Vr,-2,m ,(10) 



where the coordinates {x, y, z) fix the position of the body 
in the rotating reference frame. In equation (6), the d/d» 
means "partial derivative with respect to a coordinate x, 
y or . These derivatives (first and second) are given by 

(a, 7 e N I a -h /3 7 < 2) 



d°'x d^y d~fi 



E 

j=0 



2"+0 (n-m)! 

X C!a,l3,j Vn+a+l3+-f,m+a+l3-2j , (H) 



where 



'0 
j - kj 



(12) 



with max{0, j — a} < fc < min{/3, j}. 

Equation (4) is used to calculate the acceleration in 
the fixed reference frame. The coefficients Cn,m and Sn,rm 
which depend on the shape of the body, are provided by 
the user. 



2.2.2. Non-sphericity of the central body: MacMillan 
potential for an ellipsoid 

If the shape of the central body is approximated by a tri- 
axial ellipsoid (with semi-axes a > 6 > c, a on the x 
axis, h on the y axis and c on the z axis), the potential 
induced by this body in the rotating frame is given by 
the MacMillan formula (MacMillan 1958), reformulated by 
Garmier & Harriot (2001): 



V{x,y,z) 



where 



3n 



1 

ds 



fc2 



V/(S2 -/l2)(s2 
7,2 - &2 < q2 _ ^2 ^ 



fc2. Ai is the first elliptical 



coordinate of the point {x.y.z) calculated by taking the 
square root of the largest solution in of the equation 



The acceleration caused by this potential has been intro- 
duced in NIMASTEP (Compere et al. 2012): 



/ 



'^rot 



3 x/i y 

2 aT 7-1^ 



-6ynX 



\ 

with a = 4Ai 



i z/i Al 



/2(t)/3i/2(t) 

/■I (i + l)2 
' y_i a3/2(f)/3i/2(t) 



dt 



_1 Ql/2(t)/33/2(j) 

h'^ {t + If and /? = 4Af 



dt 



(15) 



e {t + 1)2. 

The integrals in (15) are computed with a Gauss-Legendre 
quadrature. 

To obtain the acceleration in the fixed reference frame 
we again use equation (4). 

2.2.3. Third body attraction 

The acceleration of a satellite orbiting a central point- mass 
perturbed by a third body (in a fixed reference frame linked 
to the central body) is given, for example in Montenbruck 
& Gill (2000) or Murray & Dermott (1999), by 

/ \ 



rzb = -G 



M + m 



\ Direct part 




where M, M^),, and m are the masses of the central body, 
of the third body and of the satellite. The position vec- 
tor rsi, of the third body with respect to the central body 

is given by ephemeris. These ephemeris are either highly 
accurate, like the ,JPL DE405 (Standish 1998) and INPOP 
ephemeris (Fienga et al. 2008; Manche 2010), or correspond 
to a simple Keplerian motion of the third body. 
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2.2.4. Direct solar radiation pressure 



2.2.5. Thrust acceleration 



We assume that the small body or the satellite is spherical. 

The albedo of the central body is ignored and, in a first 
step, the central body - shadowing effects are not taken 
into account. The acceleration induced by the direct solar 
radiation pressure is given by (Montenbruck & Gill 2000) 



Pq Cr 



0,0 



A 
m 



■ re 



(17) 



where Cr is the adimensionnal reflectivity coefficient that 
depends on the optical properties of the satellite surface; 
Pq = 4.56 X 10~® N/m^ is the radiation pressure for an 
object located at the distance of 1 AU; uq is a constant 
equal to the mean distance between the Sun and the 
central body, and vq is the planetocentric position of the 
Sun. Finally, the coefficient A/m is the area-to-mass ratio 
where A is the effective cross-section and m the mass of the 
satellite. This coefficient is present also in the (neglected) 
atmospheric drag effect. 

The central body - shadowing effects can be taken into 
account by multiplying formula (17) by a coefficient v. In 
this software, we offer two possibilities to calculate this co- 
efficient. In both cases we assume a spherical central body 
and we neglect its atmosphere. 

First, we assume a cylindric umbra. In this case, the Sun 
is assumed to be distant enough from the central body so 
that the solar rays are considered parallel. Two conditions 
(Montenbruck & Gill 2000) have to be respected to consider 
that the satellite is in the umbra {u = 0): 



Tq ■ r 



< and sia-ip \\r\\ < Re , 



(18) 



where ip is the angle between the vector central body to 
satellite (r) and the vector central body to Sun (r©). The 
symbol • designates the scalar product. Out of the shadow 
the coefficient v is equal to 1. 

Second, we can take into account the penumbra tran- 
sition (conical shadow). We use the formula given by 
Montenbruck & Gill (2000) to calculate the coefficient v. 
The satellite is 

— completely out of the umbra and penumbra {v = 1) if 

a + b < c; 

— completely in the umbra (i^ = 0)ifc<6— aorc<o — 6; 

— in the penumbra otherwise. 

In this last case, the coefficient f is calculated with 



where 



a arccos 



X = 



c2 + a2 - 62 



2c 



a = arcsm 



h = arcsin — - , 



arccos 



-r ■ (rQ - r) 

r\\rQ - r\\ 



(19) 

(20) 
(21) 

(22) 

(23) 



Assuming that the small body is an artificial satellite, we 
also include a constant thrust in three directions in the 
accelerations. We consider the mass of the satellite to be 
constant. The three allowed directions are the radial direc- 
tion (subscript 1), the direction of the angular momentum 
(subscript 2), and the transverse direction (subscript 3). 
The acceleration is given by (Montenbruck & Gill 2000) 



1 r X r r 



with Rq the equatorial radius of the Sun. 



r X r 

\ ^ \\r X f|| / 

where Pj is the thrust (kgm/s^) coordinate in the direction 
i, for i = 1,2,3 and the symbol x designates the vector 
product. This option has been implemented and success- 
fully used in Delsate (2011a). 

2.3. Neglected forces and effects 

At the conception of NIMASTEP, we have chosen to intro- 
duce forces by decreasing order of magnitude (for our test 
cases). Therefore some forces are not taken into account in 
NIMASTEP at the moment, but they could be inserted in 
next versions of this software. 

For example, the general relativity is neglected for two 
reasons. On the one hand, the allowed choice of the cen- 
tral body is restricted to the planets and asteroids of the 
solar system, for which the effect of the general relativity is 
weak. For instance, at the altitude (a « 12 000 km) of the 
LAGEOS satellite, the ratio between the Schwarzschild cor- 
rection and the attraction of the Earth is close to 1 x 10~®. 
This induces a secular effect on the pericenter, of this satel- 
lite, reaching 1 km after five years (Deleflie 2002). On the 
other hand, if this force is taken into account, we have to 
include some other forces at the same order of magnitude 
(see Figure 1). 

The atmospheric drag has also been neglected because, 
except for the Earth, either the models of atmosphere 
are not available or the central body has no atmosphere. 
Because our first interest was in high altitude satelhte, we 
neglected this efl^ect in the code's first conception. This is 
also the reason we did not take the tidal effect into account. 

Excluding the atmospheric drag (neglected effect) and 
the radiation pressure, the shape of an artificial satel- 
lite has an insignificant effect on its orbit (Celletti & 
Sidorenko 2008). For the radiation pressure, we considered 
a mean area-to-mass ratio. The same choice was made in 
the STELA software (see section 4.2.1): the ONES (The 
French Space Agency) decided to use an "equivalent" con- 
stant A/m. for the dynamics of non-spherical satellites. This 
is why the spacecraft is assumed to be spherical. 

Figure 1 presents graphically the order of magnitude 
(radial component) of different forces (see Capderou 2005 
and Montenbruck & Gill 2000 for simplified formulae). This 
graphic gives an idea of which forces are important or not 
with respect to the central body and the altitude of the 
satellite. This plot teaches us that to take into account the 
relativistic effect on an artificial satelhte of the Earth, we 
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Fig. 1. Order of magnitude (radial component) of the main forces acting on a small body orbiting Mercury, Earth, Mars, 
and Vesta with respect to the relative distance (ratio between the distance from the center of the central body and its 
equatorial radius) from the center of the central body. The dashed vertical lines localize the synchronous orbits. 



cannot neglected the tidal effect or, for low altitude, the 
atmospheric drag effect. For a probe orbiting Vesta, the 
relativistic effect is less than f x 10~^^ km/s^ and the effects 
of Jupiter or Sun are the strongest. 

3. Additional tools 

3.1. MEGNO: Mean exponential growth factor of nearby 
orbits 

To study the regular or chaotic behavior of the computed 
orbits, we incorporated an additional tool into NIMASTEP: 
the chaos indicator MEGNO, which stands for mean expo- 
nential growth of nearby orbits (Cincotta & Simo 2000). 
Here we introduce it briefly. 

Let us consider the dynamical system 



d , . 



flx{t)\, with X 6 



(25) 



Let (j){t) be a solution of the flow defined by equation (25), 
then the evolution of a tangent vector S{t) along (j){t) is 
given by 



S = J{<i>{t)) S{t) , with 



Jim) = ^im)- (26) 



J is the Jacobian matrix of the system (25). The MEGNO 
indicator, based on a modified time-weighted version of the 
integral form of the Lyapunov characteristic number, is de- 
fined by 



Yit) ^ - 



56 



s ds, 



(27) 



and the mean MEGNO is 

Y{i) = \ j\{s)ds. (28) 

This quantity is a characterization of the divergence rate of 
two nearby orbits. For stable periodic orbits, Y converges 
toward 0, for stable quasi-periodic orbits, Y converges to- 
ward 2 and for chaotic motion, Y tends toward infinity 
(Cincotta et al. 2003). In NIMASTEP, we adopt the same 
strategy as Gozdziewski et al. (2001) for the numerical com- 
putation of the MEGNO. 

This indicator has been successfully used in some 
papers, for instance, Breiter et al. (2005), Cincotta & 
Simo (2000), Gozdziewski et al. (2001, 2008) and Hinse 
et al. (2010, 2008). For the results concerning the MEGNO 
implemented into NIMASTEP we refer the reader to Valk 
et al. (2009a) and Compere et al. (2012). 

To integrate the system (26), we need the expression 
of the Jacobian matrix J. In NIMASTEP, the MEGNO 
indicator can be calculated for motions including the non- 
sphericity of the central body's gravity field, the third-body 
effect and the direct radiation pressure without shadow 
[v = 1). For these three contributions, the Jacobian ma- 
trix has the same structure: 







Tpart. _ dr 



I 











f 











I 







(29) 
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Because these forces are conservative, the (partial 
Jacobian) matrix JP'^'^'- is symmetric and its trace is equal 
to 0. Let us give the JP^^*- for the three cases. 

For the non-sphericity of the central body, given by the 
spherical harmonic representation, the Cunningham (1970) 
method gives 





/ d'^V 


dxdy 


d^V \ 
dxdz 




xpart. 

J J, = n 


dydx 




dydz 


, (30) 




\ dzdx 


d^V 
dzdy 


d^z / 




d^V { 


oo n 

5^ 5^ -^e {Cn,m 
71=0 m=Q 




d»d-k J 



(31) 

The second derivatives are given by (11) and this Jacobian 
matrix (30) is calculated in the rotating reference frame. 

The Jacobian matrix due to the non- uniformity of the 
central body is given by 

j-part. _ p-1 j-part. p _ pT 7-part. p ,.jr)\ 
"^pot — J^e T tig — Kg J rp tig . [ 61] 

If the central body is approximated by an ellipsoid, 
we also use equation (32) where J^*"^*' is calculated by 
(Compere et al. 2012) 



jpart. 



/ Urot _i_ f dXi f dXi f d\i 

X ^ dx dy •/ 1 9z 

V 



•'3 dx 



dy 
Jd a,, 



/2 -gj 



(33) 



with fi = fi{Xi,h,k,x) 

/2=/2(Al,/l,fc,2/) 

f3 = f3{M,h,k,z) 



3a:/i 



A2(A2-/l2)l/2(A2-fc2-)l/2 ' 



(Af - /l2)3/2(A2 _ A;2)l/2 

3zfj, 



, (34) 



(Af - /l2)l/2(A2 - fc2)3/2 • 

Xrot, ijrot and Zyot are the three components of the acceler- 
ation (15) and the partial derivatives of Ai are given in the 
appendix of Compere et al. (2012). 

The partial Jacobian matrix of the third-body effect is 
given by (Montenbruck & Gill 2000) 



jpart. 



rzbW 
1 



■ {r - rsb) (g) (r - rsbf 



■ r3b\ 



(35) 



where ® is the tensorial product and /ax a is the unitary 
matrix of dimension 3. 

Last, for the direct radiation pressure (for v = I), 
the partial Jacobian matrix is similarly calculated using 
(Montenbruck & Gill 2000) 



Jrp'^ — P@ dQ — 



dr - re) » (r- - r©)^ 



+ 



(36) 



3.2. NAFF: Numerical analysis of the fundamental 
frequencies 

Another tool directly integrated into NIMASTEP is the fre- 
quency analysis introduced by J. Laskar (Laskar 1999, 2003, 
2005). The NAFF (numerical analysis of the fundamental 
frequencies) is a very efficient numerical method for a fre- 
quency analysis, much more efficient than a classical FFT. 
Indeed, for a KAM solution, the accuracy of the frequencies 
of a signal on a time span [— T, T] is proportional to 1/T^ 
for the NAFF, using the Hanning window, while for an or- 
dinary FFT method, this accuracy is only proportional to 
1/T (Laskar 1999). 

The main purpose of the NAFF is to determine the ap- 
proximation /• {t) (with a given number N of harmonics) 
of a signal f{t) from its numerical knowledge, where both 
functions / and /* are developed in Fourier series: 

N 

f'{t) = ^ a* e*"''* is the approximation of the initial 

k=l 

oo 

signal f{t) = ^ flfc e'"''*, where fk and iv* are the real 

frequencies, when at and a' the complex amplitudes. The 
frequencies for k = 1,...,7V and their associated de- 
creasing amplitudes a* for k = 1,...,A'' are determined 

through an iterative scheme. Indeed, to determine the first 
frequency ly', one searches for the maximum of the ampli- 
tude of 

0(i.) = (/(t),exp(ii.f)), (37) 
where the scalar product {f{t),g{t)) is defined by 

{m,9it)) = ^ f_J{t)-^x{t)dt 

where g{t) is the complex conjugate of g{t) and wi. 
is a weight function, i.e., a positive function verify 



(38) 



•mg 



(39) 



Laskar (1999) recommends the following choice (Hanning) 
2P(p!^2 



x{t) 



{2p)\ 



■{l + cos{-Kt)Y . (40) 



where p is a positive integer. In practice, the algorithm is 

more efficient when p = 1 or p = 2. 

Once the first periodic term exp(i;y*i) is found, its 
complex amplitude n* is obtained by orthogonal pro- 
jection, and the process is restarted with the remainder 
fi{t) = f{t) — oJexp(ji^*i) as initial input. The algorithm 
stops when two detected frequencies are too close to each 
other, because this alters their determinations, or when the 
number of detected terms reaches a maximum set by the 
user. When the difference between two frequencies is larger 
than twice the frequency associated with the length of the 
total time interval, the determination of each fundamental 
frequency is not perturbed by the other ones. 

To our knowledge, the first uses of NAFF in the 

geodesy context were published by Delsate & Lemaitre 
(2009), Delsate et al. (2008, 2009) and Valk et al. (2009a). 
In Valk et al, (2009a), the authors (using NIMASTEP and 
NAFF) computed the frequency of the resonant angle of 
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the geostationary orbit. Their results agree weU with the 
theory of Valk et al. (2009b). 

It is possible to use NAFF for other purposes than the 
decomposition of a signal into frequencies: it can be used 
as an indicator of dissipation of the main frequencies and 

therefore as a chaos indicator. This can be achieved in two 
ways. The first possibility is to calculate the variations of 
the main frequency (let us call it v) of a signal with respect 
to the initial conditions (here e.g. x and y). We use the 
frequency numerical second derivative, 6Sv, given by the 
formula (Laskar 1993) 

56i'{x,y) = [^{x.y) — 2v{x — Ax,y) + u(x — 2Ax,y)\ 

+ \i^{x,y) - 2iy{x,y - Ay) + v{x.y - 2Ay)|. (41) 

With \og(56i'{x,y)), we can detect the irregularities in a 
frequency map. The second possibility is to calculate the 
variations of the main frequency of a signal with respect 
to the time; we compute a first value v^^' of the main fre- 
quency of the signal on a time span Ti (e.g. from —T to 0), 
then we compute a second value f^^' of the same signal on a 
time span (e.g. from to T). Consequently, we estimate 
the diffusion rate with respect to the time by 

Because the two indicators (41 and 42) give approximately 

the same results, they can be used interchangeably (Laskar 
1999). The computation of the second-order derivatives 
usually requires very small stepsizes Ax and Ay, while the 
computation of the time diffusion requires more iterations. 

These two methods were successfully used with 
NIMASTEP in Lemaitre et al. (2009) and Delsate et al. 
(2010). 

4. Validation 

We performed very many tests to prove that the algo- 
rithms programed in NIMASTEP were correctly imple- 
mented. Because it is not possible to describe all tests, we 
chose to show some validations of the code related to ap- 
plications (already published or not). We present different 
comparisons with an external full numerical integration, 
semi- analytical studies, and some known analytic solutions. 

We usually do not have an easy access to real data, 
therefore comparisons with these data are difficult to per- 
form. In addition, it is much easier and rigorous to compare 
the performance or results with that of other codes or with 
analytical results. The main reason is that for these data 
we have all informations on the approximations. Real data 
are not ideal to validate a code based on an approximate 
model of reality. 

4.1. Comparisons with full numerical integrations 

For a direct comparison (section 4.1.1), we compare 
NIMASTEP with "real" data and with a software written 
by A. Rossi called R-ISTIs in the following. For an indirect 

comparison (section 4.1.2), we compare measures (as the 
evolution of the MEGNO or the semi-major axis variation) 
computed through integrations obtained by NIMASTEP 
with external published results. 



4.1.1. Direct comparison 

The first validation is the comparison of an integration per- 
formed by NIMASTEP with "real" data and with the re- 
sults of R-ISTIs [personal communication of A. Rossi) us- 
ing the same initial conditions, parameters and, forces. 

First, we briefly describe the TLE (two-line elements) 
data base (because it is used by R-ISTIs and to compare our 
software with "real" data.) A two-line element^ is a data 
format used by NASA, to describe the orbital elements of 
the orbit of an Earth satellite. The elements in the TLE sets 
are mean elements calculated to fit a set of observations us- 
ing a specific model: the SGP4/SDP4 orbital model (Hoots 
& Roehrich 1980). The mean orbital elements produced by 
the SGP4/SDP4 orbital model are Earth-centered-inertial 
coordinates with respect to the true equator and the mean 
equinox of epoch. They do not include the effect of nuta- 
tion. 

It is important to note that the arithmetic or geometric 
means of osculating data do not have the same value as the 
TLE. In the same way, different orbital models give different 
mean elements. Indeed, the TLE of a satellite are neither 
the real osculating Keplerian elements, nor the mean ones. 
However, we consider that the TLE gives a good idea of the 
real orbit. 

Second, we present a short description of the method 
used in R-ISTIs. R-ISTIs converts the TLE of the satel- 
lite Etalon-1 (catalog number 19751) to osculating elements 
using the SGP4 theory (Hoots & Roehrich 1980). Then, 
the integrations are realized thanks to a special pertur- 
bation propagator with an integrator based on Cowell's 
method (Lyddane & Cohen 1962; Yoon et al. 1997) for 
the numerical integration of the equations of motion. The 
force model includes the Earth's gravity potential, the luni- 
solar gravitational perturbation, the solar radiation pres- 
sure with eclipses, and several atmospheric density mod- 
els for air drag computation. This propagator was assem- 
bled in Pisa, at ISTI^, based on the NASA program called 
ASAP-artificial satellite analysis program (Kwok 1987) and 
used, for example, for some of the propagations in Rossi 
(2008). At the .Julian date 2448135.5, the initial condi- 
tions of this simulation, calculated by A. Rossi (R-ISTIs), 
are a = 25 501.226477 km, e = 0.642 773427 x lO'^, 
i = 1.132 591133 rad, = 2.726 705 844 rad, w = 
4.284489314 rad and M = 0.243368 293 rad. 

The model of forces, used in NIMASTEP for this test, 
includes the Earth's attraction developed in a spherical har- 
monic expansion up to degree and order 20 (the Earth grav- 
ity field adopted is the EGM96 model, Lemoine 1987), as 
well as the combined attractions of the Sun and the Moon 
(ephemeris DE405 given by the JPL, Standish 1998). The 
perturbing effects of direct solar radiation pressure are also 
taken into account for an area-to-mass ratio fixed at lO"'* 
m^/kg with conical shadow. At this altitude, the atmo- 
spheric drag is insignificant. 

We chose the ABMIO integrator with a step size of 
200 seconds of time. Figure 2 (upper panels) shows the evo- 
lution of the Keplerian elements for both softwares showed 
in the (sinz/2cosn,sinj/2sini2) and (ecosw, esinw) 
planes. We also overlaid the data firom the TLE data base 

^ The TLE are freely available on https : / /www . space-track . 
org or http : // celestrak . com 

^ Istituto di Scienza e Tecnologie dell'Informazione "A. 
Faedo". 
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Fig. 2. Comparison of NIMASTEP witli R-ISTIs and witli tiie TLE data of Etalon-1. Two upper panels: projection of 
the orbits performed witii NIMASTEP (light blue dots), R-ISTIs (red circles), and TLE (blue stars) on the phase spaces 
(sini/2cosfl,sini/2sinf2) and (e cose*;, e sincj). We added 2.7 x 10~^ to esinw of the TLE to obtain a good match. Six 
lower panels: differences between NIMASTEP and R-ISTIs orbits. 



with a 2.7 x 10 * correction added to esintj to obtain a 
good agreement. This correction is necessary because the 
initial conditions taken by R-ISTIs and NIMASTEP are 
osculating elements and are not exactly the same as the 
TLE, probably deriving from the definition of the TLE. 
Remember that TLE are sets of mean elements but not 
equal to the mean of the osculating Keplerian elements. 
This is also why we do not present any differences (which 
will be irrelevant) between TLE and NIMASTEP results. 
Both softwares give results that agree with the TLE except 



for the eccentricity and the argument of pericenter, where 
the TLE data were parallelly shifted up. 

There are some differences between the orbit of R- 
ISTIs and NIMASTEP, but the errors (six lower panels 
of Figure 2) are very low, which proves the reliability of 
NIMASTEP in this context. The amplitude of the errors 
increases with respect to the time, but the mean error stays 
quasi-constant except for the mean anomaly. This can be 
explained by a difference of formulation of the Earth rota- 
tion between both softwares. 
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4.1.2. Indirect comparison 

Tricarico & Sykes (2010) and Delsate (2011a) studied the 
main gravitational resonances (resonances between the rev- 
olution ol the satellite and the rotation ol the asteroid, 
hereafter denoted by GR) appearing around the asteroid 
(4) Vesta. A part of these works consists in calculating the 
probability of capture of the Dawn spacecraft mission in 
the 1:1 GR of Vesta. 

Tricarico & Sykes (2010) performed this test through a 
full numerical software. They simulated a slow descent of 
Dawn from an initial radius of 1 000 km using a 20 mN 
of thrust and obtained a probability of capture equal to 
1/12 = 8.33%. They also observed that higher thrust- 
ing (until 50 mN) can exhibit a similar behavior. Taking 
the same initial conditions (semi-major axis) and forces 
model, Delsate (2011a) performed numerical integrations 
with 50 000 random mean anomaly with NIMASTEP for 32 
different thrusts from 20 mN to 36 mN. He obtained a mean 
probability of 8.262%, which agrees well (indirectly) with 
the results obtained by both softwares. That also shows 
that NIMASTEP is well adapted to this type of study. 

4.2. Comparisons with semi-analytical theories 

As a second validation of our software, we applied 
NIMASTEP to the typical case of an high-altitude aban- 
doned space debris and compared otir results with two 
semi-analytical theories: the STELA software and a semi- 
analytical study of Valk et al. (2008, 2009b). 

4.2.1. STELA software 

The semi-analytic tool for end of life analysis (STELA)"* is a 
semi-analytical software designed by the ONES, which com- 
puted efficient long-term propagations of LEO (low earth 
orbit), GEO (geostationary earth orbit), and GTO (geo- 
stationary transfer orbit) orbits based on semi-analytical 
models. It allows the computation of "graveyard orbits" 
and tests the safety of protected zones. 

STELA computes the mean orbit parameters by a semi- 
analytical model (without short periodic motion) at each 
integration time step. To compute the osculating parame- 
ters at a given time, STELA uses the mean (averaged) el- 
ements and artificially adds the short periodic effects. The 
short periodic computation model contains the influence 
of the oblateness (C20) of the Earth, the solar and lunar 
gravity, and the solar radiation pressure. The mean motion 
(middle- and long-term evolution of the orbital parameters) 
includes a complete degree and order 4 model of the grav- 
ity field for the Earth, the solar and lunar gravity and the 
solar radiation pressure. The Earth's shadow is cylindric 
and the object is assumed to follow a quasi-circular orbit. 
The integrator is a sixth-order Runge-Kutta method with 
a step size of one day. More explanations are available, for 
instance, in Deleflie et al. (2005, 2010), Fraysse (2011), and 
Le Fevre et al. (2011). 

We compared the osculating motion given by 
NIMASTEP (with the same contributions) with the re- 
sult of STELA. We chose the following initial conditions: 
a = 42164.140 km, e = 0.001, i = 0.1 rad, and Q = lj = 
M = rad. The initial time at epoch is 2451350.5 JD and 

http : / /logiciels . cnes . f r/STELA/ f r/logiciel . htm 



the area-to-mass ratio is equal to 0.1 m^/kg with a reflec- 
tivity coefficient equal to 1. The integrator is ABMIO with 
a step size of 864 sec. Figure 3 shows the difl'erences be- 
tween both integrations (in absolute error) . We notice that 
after 100 years, the maximal error is 1.911 km for the semi- 
major axis, 2 X 10~^ for the eccentricity, 30.94 arcmin for 
the inclination, and 214.51 arcmin for the ascending node. 
The amplitude of the error on the argument of pericenter 
may be caused by a small phase difference in this motion. 
Another possible explanation is that the method used by 
both softwares to obtain the pericenter from the cartesian 
coordinates or from the equinoctial elements are different. 
Moreover, the eccentricity is very low, inducing an eventual 
difficulty to determinate the pericenter angle. 

The errors are very low for this long time of integration 
and we conclude that the results NIMASTEP agree well 
with those of the STELA software. 



4.2.2. Vallc semi-analytical study 

The second validation by a semi-analytical study is 
achieved by comparing our results with the mean motion 
given by the theory of S. Valk, who studied the geostation- 
ary Earth orbits (Valk et al. 2008, 2009b). 

Directly comparing the results of the osculating numer- 
ical integration (NIMASTEP) with a semi-analytical solu- 
tion is not entirely consistent. Indeed, the osculating el- 
ements and the mean ones are not directly comparable. 
However, we consider that our main purpose is to give a 
first order comparison. 

We start with a short description of the semi-analytical 
theory; for more details we refer the reader to Valk et al. 
(2008, 2009b). The framework is a mean motion semi- 
analytical theory: only the long-term and secular effects 
are derived. In other words, the resulting theory does not 
include any short-term effects. In practice, the theory con- 
sists of the numerical integration of the filtered equations 
of motion over the short periods. 

First, the approach is to choose a nonsingular and 
canonical set of variables, namely the Poincare variables: 



xi = v /2 (L - G) sm{-Lj - f2) , 
X2 = v /2 {G - H) sin(-r2) , 
= v^2 (G - H) cos(-n) , 
a;3 = A = 9. + u) + M , 
Xi — -sjl [L — G) cos(— w — f2) , 
x^ = L, 



where L = yfjm, G = L\j\ — and H = Gcosi are the 
Delaunay elements. Valk et al. (2008, 2009b) expanded the 
Earth potential in powers of the eccentricity and of the 
inclination, truncateed the development at a low order to 
obtain an Hamiltonian formulation. Some other perturba- 
tions were added, such as the direct radiation pressure (Valk 
et al. 2008) and the solar and lunar gravitational effects. 
Second, the final expansion of the Hamiltonian is set in 
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Last, they averaged the disturbing function over the fast 
variable, namely the sidereal time 0, to the first order by 
dropping the fast periodic terms in the trigonometric ex- 
pansions. The Hamiltonian is written as a Poisson series 
using the symbolic software MSNAM, the series manipula- 
tor of Namur (Henrard 1986). 



For our validation, let us consider the dynamical evo- 
lution of a theoretical high-altitude space debris. The per- 
turbations taken into account are the oblateness (C20) of 
the Earth and the third-body perturbation induced by the 
Moon. The integration time is equal to 100 years. The 
entry-level step size used in the osculating numerical in- 
tegration with NIMASTEP is fixed at 200 seconds (with 
the RK4 method), whereas Valk et al. (2008) define a one 
day step size to integrate the averaged system of equations 
(also with the RK4 method). 

To quantitatively estimate the accuracy of our soft- 
ware, Figure 4 shows the differences between both orbits. 
First, the differences remain small although the chosen time 
of integration is very sizable. Indeed, the computed root 
mean square (RMS) can be considered to be very small 
because the comparison is made without any preliminary 
fitting of the mean initial conditions (Deprit 1969; Henrard 
1970: Valk et al. 2009a). Moreover, the RMS on both the 
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44 164 km, eo = 0, iq = t<Jo = -^0 I'ad. Initial time at epoch is 25 January 1991. 



argument of perigee and the longitude of the ascending 
node are mostly influenced by the singular transformation 
when projecting the integrated non-singular state vector 
into Keplerian orbital elements. 

4.3. Comparison with results of analytical studies 

To validate our software yet again, we compared some 
results of NIMASTEP with predictions of analytical stud- 
ies. Some previous articles have shown a good agreement 
between analytical results and numerical integrations 
provided by NIMASTEP. In this paper, we only give the 
context of each analytical study and how our software is 
validated. The reader is refered to these articles for more 
details on the results and comparisons. 

Valk et al. (2009a,b) and Lemaitre et al. (2009) studied 
the motion of space debris with or without high area-to- 
mass ratio close to the geostationary Earth orbit (GEO). 
They take into account the ellipticity of the equator (C22) 
and the direct solar radiation pressure. At first, Valk et al. 



(2009a) assumed A/m = and localized numerically (with 
NIMASTEP) the position of the GEO equilibrium. The 
analytical study of Valk et al. (2009b) yielded the same lo- 
cation. Second, Valk et al. (2009a) performed a frequency 
analysis (with NAFF in NIMASTEP) of the resonant an- 
gle evolution. They obtained a fundamental frequency of 
libration differing from 0.44% with the analytical estima- 
tion (Chao 2005; Valk et al. 2009b; Vallado 2001). Third, 
increasing the area-to-mass ratio, Valk et al. (2009a) dis- 
covered (always helped by NIMASTEP) a web of secondary 
resonances close to the GEO resonance. Delsate & Lemaitre 
(2009) and Lemaitre et al. (2009) analytically explained, lo- 
cated and listed some of the fundamental frequencies. 

Delsate et al. (2010) studied the stability of a massless 
probe orbiting an oblate central body perturbed by a third 
body. They analytically determined the location of frozen 
orbits (e = w = 0). They also computed the periods of 
the equilibria. Moreover, for the particular case of a probe 
around Mercury, they used NIMASTEP to validate their 
results. They satisfactorily recovered the locations, orbits 
and periods of the equilibria numerically. 
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As previously said (section 4.1.2), Delsate (2011a) stud- 
ied the main GR appearing around the asteroid (4) Vesta. 
NIMASTEP was used to obtain a map of the radius vari- 
ation of a spacecraft bringing the GR to the fore. These 
numerical results agree well with those of Tricarico & Sykes 
(2010) and with the analytical study. 

Compere et al. (2012) studied the evolution of a small 
body in rotation around a rigid ellipsoidal asteroid in con- 
stant rotation. Simulations were made with NIMASTEP 
using the MacMillan (1958) potential for the influence of 
the shape of the asteroid on the small body and the chaos 
indicator MEGNO. The authors presented stabihty maps 
showing extremely stable conic-like curves corresponding to 
GR. An analytical model, based on a simplification of the 
MacMillan potential, was created to validate the numeri- 
cal results. Here again, the analytical model confirms the 
results performed with NIMASTEP. 



5. Summary and perspectives 

The dedicated software NIMASTEP allows one to numer- 
ically integrate the osculating motion of an object (artifi- 
cial or natural satellite, space debris, small natural body) 
orbiting a telluric planet, a (dwarf-)planet, or an aster- 
oid of the solar system. The equations of motion, in the 
Newtonian formalism, can include different types of pertur- 
bations, such as the influence of one (or more) third body 
(bodies), of the non-spherical shape of the central body, 
of the radiation pressure, or even of a thrust acceleration. 
This software includes two supplementary tools to analyze 
the orbits: the chaos indicator MEGNO and the frequency 
analysis NAFF. 

It is flexible and modular, which allows very many op- 
tions (forces, central bodies, object characteristics) to be se- 
lected or disabled easily and, in the future, to add some pos- 
sibilities and the neglected forces. By this extensibility, this 
software could be used by many members of the dynami- 
cal or geodesian community. It has been successfully val- 
idated by several astrophysical interest tests, comparisons 
with other numerical integrations, with semi-analytical re- 
sults, and finally with analytical studies. 

The code is not available yet on the internet, but 
a web site provides information for contacting its au- 
thors and fo obtaining (in the near future) the software: 
http : //www . f undp. ac . be/en/research/projects/ 
page_view/10278201/. 

The most important future additions scheduled for 
NIMASTEP are (I) to include the atmospheric drag of the 
central body, (II) to simultaneously calculate the motion of 
more than one (interacting) satellites, (III) to include the 
general relativity (Will 2006), (IV) to include tidal effects, 
(V) to include the attitude motion of the satellite, (VI) to 
modelize the gravity field by other methods, and (VII) to 
add YORP effects. 
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